Direct evaluation of large-deviation functions 
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We introduce a numerical procedure to evaluate directly the probabilities of large deviations of 
physical quantities, such as current or density, that are local in time. The large-deviation functions 
are given in terms of the typical properties of a modified dynamics, and since they no longer 
involve rare events, can be evaluated efficiently and over a wider ranges of values. We illustrate the 
method with the current fluctuations of the Totally Asymmetric Exclusion Process and with the 
work distribution of a driven Lorentz gas. 
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In the last few years there has been a renewed in- 
terest in the theory of large deviations of nonequilib- 
rium systems, with the development of general results 
concerning the fluctuations of soft modes Q, of non- 
trivial and rich analytic solutions of explicit models 
HSSEEHSHOjii and with the discovery 
of strikingly simple an d g eneral nonequilibrium rela- 
tions [H tH H Q E El [H El [3 (the Fluctuation 
Theorem, Jarzynski's relation) obeyed by work fluctua- 
tions. Perhaps for the first time, we are gathering a few 
glimpses of truly general features of macroscopic systems 
well out of equilibrium. 

The large-deviation function plays an essential role in 
the investigation of nonequilibrium systems — a role akin 
to the free energy in equilibrium ones. When available 
techniques do not allow for an exact evaluation of this 
function, one turns to simulations: but direct numerical 
simulation of large deviations is hard, since, by defini- 
tion, they are rare. In equilibrium, this difficulty is often 
overcome by introducing biased (non-Boltzmann) sam- 
pling |2jj. Here we show that a similar strategy can be 
introduced in systems out of equilibrium, in order to eval- 
uate the large deviations function for quantities that are 
local in time, although not necessarily in space. We find 
that, in nonequilibrium, it is necessary not only to bias 
suitably the dynamics of the system, but also to intro- 
duce a process by which images (clones) of the system 
reproduce or die, a technique inspired by the Diffusion 
Monte Carlo method [U of quantum mechanics. In the 
present work, after deriving the general formalism, we 
show its effectiveness by applying it to two nonequilib- 
rium processes: a stochastic one, the Totally Asymmetric 
Exclusion Process (TASEP), and a deterministic one, a 
driven Lorentz gas. Our algorithm allows to compute 
the probability of obtaining a temporary large deviation 
(compared to the typical) value of the current in the first 
example and of the dissipated work in the second one. 

We consider the general setup of a system evolving ac- 
cording to Markovian dynamics. Let C,C denote two 
configurations in the and let Uc'c be the transition ma- 



trix of the discrete (eventually continuous) time dynam- 
ics. Denoting by Pc(t) the probability of being in the 
configuration C at time t, one has 



P c >(t + l) = J2 u c'cPc(t). 



(1) 



In a time interval of length T, a path Co, C\, . . . , Cr in the 
configuration space, starting from a fixed Co, will have the 
probability 



Prob[C , Ci, . . . ,C T ] = U Ct 



Ct- 



Uc 2Cl -U Cl c . (2) 



ditive in time, i.e., which can be written as Qt — 
Y^t=o J{Ct+i,Ct). For example, for a transition C — > C 



We shall consider physical quantities Qt that are ad- 
ditive in time, i.e 

X)t=o J{Ct+i,C t ). 
in a lattice system: 

Jcc — 



1, a particle jumps to the right; 
0, nothing happens; (3) 
— 1, a particle jumps to the left. 



We are interested in calculating the probability of having 
a current Qt in the time interval T, i.e., a current q = 
Qt/T per unit time. Denoting with angular brackets the 
average over trajectories, we have 



JCtCt- 



+200 



2-ki 



(4) 



where we used the integral representation of delta func- 
tion and we have defined 



7>(A) _ /MJc 1 c +-+Jc T c T _ 1 ) 



C U ...,C T 



Udc e 



(5) 

M ^CxCqH VJc T c T _ 1 ) 



In the limit T — * oo, by applying the steepest de- 
scent method to Eq. (4) (and assuming that the imag- 
inary contour line can be deformed to the real line), 
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one obtains that the large deviation function f(q) = 
liniT^oo ]&P(Qt)/T and /z(A) are Legendre transforms 
of each other: 



f(q) = max[M(A) 



(6) 



so that q = //(A*) where A* is the saddle. We introduce 
the bias Q of the original measure by defining the new 
matrix 



Ucc = e AJc 



so that 



E 

d,...,e T 



CtCt- 



'Ucc, 



Udc - 



(7) 



E 

Ct 



u 1 



CtCq 



(8) 

Introducing the spectral decomposition U = 
J2j e Aj | A^) {A^ | where we assumed a complete biorthog- 
onal set of the matrix U exists, i.e. U\A R ) = e Aj \A R ) 
{A-j\U — e Aj (A-j\ and denoting by e A the eigenvalue 
of U with the largest real part, and by \A R ), \A L ) the 
corresponding right and left eigenvectors, we have, for 
large times T, 



(C T \A R )(A L \C ) e 



TA 



(9) 



so that A = /i(A). In order to compute n(X), one possibil- 
ity is to perform path-sampling over the trajectories with 
weight J7J). Such a procedure has been proposed in the 
context of the work distributions |22T | on nonequilibrium 
trajectories. In this paper we propose a different strat- 
egy: the idea is to define a new effective dynamics whose 
expectation values directly give the large deviations. As 
we shall see, the new dynamics involves the parallel evo- 
lution of clones which reproduce and die, a procedure 
inspired by the "Diffusion Monte Carlo" method of sim- 
ulation of the Schrodinger equation |2*lj . In order to write 
eq. JHJ as expectation on the new dynamics, let us put 
Kc = X)c Ucc, an d define the stochastic matrix 



U c , c = UacK c \ 
We now have, instead of eq. ©, 



(10) 



'"■■■'■= > u CtCti K Ct _ 1 



■■U'c iCq K Co . (11) 



This can be realized by considering an ensemble of L 
copies ( "clones" ) of the system, and by successively go- 
ing, for all of them, through a process defined by the 
following three steps: 



A cloning step: 



P c (t + 1/2) = K c Pc(t), 



(12) 



where the configuration C of the selected copy gives 
rise to G identical clones, G = [Kc] + 1 with prob- 
ability Kq — [Kc], and G — [Kc] otherwise ([x] 
denotes the integer part of x). If [Kc] = 0, the 
copy may be killed and leave no offspring. 

• A shift step without cloning of all the offspring of 
C with the modified dynamics V 



(13) 



• An overall cloning step with an adjustable rate 
Mt = L/(L + G) (at each time the same for all 
configurations), so as to keep the total number of 
clones constant. This amounts to multiplying U by 
Mt times an identity, at each time. 

It is easy to see that, in the long-time limit, the com- 
pensatory factor gives us /z(A) through 



- ln[M T • • • M 2 ■ Mi] = 7>(A). 



(14) 



Remark 1: We note that, if the quantity Qt, whose 
deviations we wish to compute, depends on a single con- 
figuration rather than on a pair of configurations, such 
as for the density Qt = Y^t=i p{Ct), the same derivation 
goes through with the substitution Jcc — * p(C). 

Remark 2: The configurations obtained in the course 
of the simulation are representative of the typical ones 
at the end (t = T), rather than within (0 < t < T) the 
interval of time T during which the large deviations are 
observed. (Their probabilities are proportional to(C\A R ) 
and (A L \C)(C\A R ), respectively). 

We now turn to two examples: the Totally Asymmetric 
Exclusion Process, and the Lorenz gas. 

A stochastic system: The Totally Asymmetric Exclu- 
sion Process (TASEP). The TASEP IQj consists of par- 
ticles on a ring with discrete sites with occupancy zero or 
one. A given particle chosen at random does not attempt 
to move with probability (1 — a), and with probability 
a attempts to move to the right and succeeds if the cor- 
responding site is empty. The parameter a can be made 
small to approach the continuous time limit. Here we 
shall set it to unity. Let us denote by Xq the number of 
different configurations that can be reached by making 
a one-particle move (1PM) from C. Then the non-zero 
entries of Ucc are given by 

TT ( a/N, if C -> C is a 1PM; . 

Uc ' c -\l-(X c a/N), \iC'=C. (15) 



This implies for U 

' ae x /N, 



Ucc 



if C C is a 1PM; 



1 - (X c a/N) if C =C. 



(16) 



Thus, for a configuration C with Xq mobile particles, we 
have 



^c = l + -^(e A -l), 



(17) 
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and finally 

, ( (ae x /N)/K e , if C -> C is a 1PM, 

^ c ' c " \ (1 - X c ,a/N)/K Cl if C = C. 

(18) 

Thus, with probability (1 — Xca/N)/ Kq no move is 
made; otherwise we move a particle randomly chosen 
with uniform probability among the Xq mobile particles. 




FIG. 1: Space-time diagram for a ring of N = 100 sites. Top: 
A = —50 and density 0.5; the shock is dense and does not 
advance. Bottom: A = —30 and density 0.3; the shock drifts 
to the right. 

In Figure H (top) we show a space-time diagram of the 
system with N = 100 particles, density 0.5 and A = —50. 
The simulation was done with L = 1000 clones, each of 
them initialized with random (uniform) occupancy num- 
bers, such that the configuration has density 0.5. We 
notice that the configurations rapidly become inhomoge- 
neous, exhibiting an alternation of a regions with high 
density with regions of slow density, as in traffic jams or 
in shock waves. The high-density regions eventually coa- 
lesce into a single one. The figure does not quite represent 
the evolution of a shock (because of Remark 2 above), but 
rather the configuration at the end of the time-interval 
for time intervals ending at progressively longer times. 
As predicted by the theory for this value of the density, 
the shock does not drift, although different initial condi- 
tions lead to different shock positions. Bottom of Figure 
n shows the case A = —30, and density 0.3: we see that 
the shock has a net drift to the right, again as predicted 
by the theory [3j . Finally, in figure [21 we show the nu- 
merical results obtained for /x(A), and compare them to 



the analytic ones of Ref. [3j. The agreement is excellent, 
and the numerical effort corresponds to tens of minutes 
of a pe 




X 



FIG. 2: Plot of ju(A) vs. A for the TASEP at density one-half. 
Numerical results and analytic results of ref. |3J, with points 
and full line, respectively. 

A deterministic system: The Lorentz Gas and the 
Gallavotti- Cohen theorem. This system consists of a 
number of particles (in our case only one) moving inside 
a billiard as in figure |3 with periodic boundary condi- 
tions. The particle is under the action of a force field E, 
and is subject to a deterministic thermostat that keeps 
the velocity modulus constant \v\ = 1. Between bounces, 
the equations of motion are: 

Xi = -Ei + j(t)±i, z = l,2; 
7 (f) = ^Ekxi. (19) 

i 

We wish to compute the generating function of the dis- 
sipated work Qt — Jq l(t) dt, and check the Gallavotti- 
Cohen theorem, which states that P(Qt) / 'P(-Qt) — 
cxp(Qt), which is equivalent, thanks to Eq. to the 
symmetry of around A = — \. 

The dynamics is deterministic, and hence cloned sys- 
tems will evolve together and perform a poor sampling. 




FIG. 3: The billiard. The radii are Ri = 0.39, R 2 = 0.79. 
We also show an example of trajectory for the external field 
£=(1,0). 
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FIG. 4: Plot of p(A) vs. A for the driven Lorentz gas. Data 
for E = (E, 0), E = 1,2 and noise intensity A = 10~ 3 , 1CT 4 . 
The Gallavotti-Cohen theorem implies the symmetry around 
A = —1/2. The best fit (continuous lines) is quadratic for 
E — 1 (gaussian behavior), and a 4-th order polynomial for 
E = 2. 



To get around this problem, we introduce a small stochas- 
tic noise, and check the stability of results in the limit 
of small noise. We evolve the system for macroscopic in- 
tervals T, and clone with a factor K t — e XCT ^ t \ where 

t-\-T 

IV(t) = J t j(t) dt is the total dissipated work over 
the interval. Before each deterministic step of time T, 
clones are given random kicks of variance A in position 
and/or velocity direction. The time-interval T and the 
noise intensity A are chosen so that twin clones have a 
chance to separate during time T , and this depends on 
the chaotic properties of the system. In the present case, 
we checked that 0.1 < T < 1 allows for a few collisions, 
which guarantees clone diversity for 1CP 3 < A < 10 . 

In Fig. 01 we show the results of //(A) for —3 < A < 2, 
and for E = (E, 0) with E = 1, 2, corresponding to very 
large current deviations. 

Diversity in a population of reproducing units is main- 
tained by a balance between the natural loss due to sam- 
pling fluctuations and the increase introduced by mu- 
tations, represented in our case by noise [23}. Thus, if 
the noise level is too small in the billiard case, most of 
the clones correspond to too close configurations, and 
our results become noisy and unreliable. The same phe- 
nomenon explains why all clones exhibit shocks in essen- 
tially the same position for any given run in the TASEP 
(since they share a common ancestor) : but we found that 
in this case the phenomenon poses no problem for the 
sampling, since the current does not depend on the posi- 
tion of the shock. 

In conclusion, we have shown that sampling methods 
based on a modified dynamics with clones can be used 
to efficiently compute the large deviations function, in 



times and within ranges of values that cannot be reached 
in a direct simulation. 
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